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THERMOELASTIC DAMPING AND ITS EFFECT ON FLUTTER OF 
STRESSED PANELS SITUATED IN A SUPERSONIC AIRFLOW 

By R. C. Shieh* 

Langley Research Center 

SUMMARY 

The effects of material damping on flutter of stressed rectangular panels are studied 
within the context of linear thermoelasticity theory. The closed-form expression for the 
thermoelastic (material) damping coefficient is obtained as a function of frequency, panel 
temperature and dimensions, and material properties. The solution of the stability 
boundary-value problem is obtained by use of a generalized Galerkin method in the cross - 
stream direction which reduces the governing partial differential equations to a system 
of ordinary differential equations in the streamwise direction. These equations are then 
solved exactly. Numerical results are given for the thermoelastic damping coefficients 
and for the flutter speeds of partially and fully clamped panels subjected to midplane 
stress. 


INTRODUCTION 

The effects of in-plane stress on panel flutter have been actively studied by many 
authors in recent years. Extensive references on the subject may be found in review 
articles (for example, refs. 1 to 5). The agreement between experimental and theoretical 
results has often been poor because of the difficulty in defining or measuring certain 
parameters (such as material and support damping, in-plane stresses, boundary-layer 
thickness, and support conditions) affecting the flutter boundary of a given panel. Inclu- 
sion of structural damping in theoretical flutter analyses of stressed panels is important 
since it may improve greatly the agreement between theoretical and experimental results. 
In particular, the paradoxical phenomenon of zero flutter speed predicted by theory which 
neglects structural damping can be removed (see refs. 6 and 7). As shown in reference 8, 
this paradoxical phenomenon can also be removed by accounting for the effect of initial 
imperfections of the panels. However, this approach is less desirable from a computa- 
tional standpoint since in addition to the difficulty of defining "initial imperfections," non- 
linear theory is required for the analysis. 
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structural damping for panels consists of both material damping and frictional 
damping at the panel supports. At present, the most widely used material -damping 
models in flutter analyses are either the linear viscous or hysteresis types; support 
damping has not been considered. Experimental evidence shows, however, that material- 
damping models of these types are not realistic for a wide class of metals (see ref. 9) 
including aluminum (see ref. 10). However, material damping can be closely predicted, 
based on the theory of thermoelasticity by including thermomechanical coupling, as indi- 
cated by the good correlation between experimental and theoretical results for material 
damping for beams (see ref. 11). Since structural damping has a significant effect on 
panel flutter, it seems worthwhile to study the damping mechanisms for panels and the 
effects of this damping in a more rigorous manner than in previous works through the use 
of thermoelasticity theory and through a more rational approach to the application of 
material and support damping. In the present paper the thermoelastic material damping 
mechanism, as well as the effects of such damping on panel flutter, is considered as a 
first step toward the aforementioned goal. Since thermoelastic damping is frequency 
dependent (ref. 11), inclusion of only this type of damping in panel flutter analyses cannot 
completely remove the contradictory phenomenon of zero flutter speed predicted by theory 
neglecting damping for certain special cases (ref. 6). In order to eliminate such a con- 
tradictory phenomenon completely, either frequency -independent structural damping or 
initial imperfections (see ref. 8) must be included in the theoretical analysis. 

In the present analysis a fully or partially clamped (including simply supported) 
thin rectangular panel (plate) is assumed to be in a state of equilibrium under the action 
of uniform in-plane stresses around all edges of the panel. All surfaces of the panel are 
assumed to be thermally insulated for calculations of the thermoelastic damping param- 
eters since, as shown in reference 11, the correlation between experimental and theoret- 
ical damping results based on such a thermal boundary condition is definitely better than 
that of an isothermal one. The lateral aerodynamic force induced by the supersonic air- 
flow is assumed to be given by the two-dimensional quasi-steady aerodynamic theory. 

The formulation of the coupled thermoelastic boundary -value problem is first presented. 

A closed-form expression for the thermoelastic material damping coefficient is then 
obtained, and the coupled governing equations are thus uncoupled. The uncoupled equa- 
tions are further reduced to a system of ordinary differential equations by the use of a 
generalized Galerkin technique in the cross-stream direction and the resulting equations 
are then solved exactly. Numerical results for the thermoelastic damping coefficient 
and for critical flutter speeds, together with a discussion of the results, are given. 
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SYMBOLS 

The units used for the physical quantities of this paper are given both in the Inter- 
national System of Units (SI) and in the U.S. Customary Units. Factors relating the two 
systems are given in reference 12 and those used in the present investigation are pre- 
sented in appendix A. The measurements and calculations were made in U.S. Customary 
Units. 

a panel length 


^mn>*^mn 

•^mnj^mn 


matrices defined by equations after equation (BIO) 


panel width 


speed of sound in undisturbed air 


specific heat at constant deformation of elastic solid 


plate flexural modulus, D = 


i(l - v^) 


Young's modulus 


lateral aerodynamic force (see eq. (4)) 


thermomechanical coupling function, g = g{p,v,\T) 


aerodynamic damping coefficient 


thermoelastic damping coefficient for a beam. 




thermoelastic damping coefficient (imaginary part of g) 


panel thickness 


function defined by equation (B4) 


3 


Kx,Ky 


k 

kx,ky 

M 


M»p 

N 


Nx,Ny 


P 


Pb 

Px’Py 


q 

R 

T 

Ta 

To 

Ts 

t 


integers 

rotational spring constants 
coefficient of heat conduction 

nondimensional rotational spring constants; kjj = K^sl^D, ky = Ky^D 
Mach number of undisturbed airflow 
thermal moment (see eq. (1)) 

number of terms of trial function in cross -stream direction (see eq. (B5)) 
in -plane thermal load defined by equation (8) 

in-plane compressive forces per unit length of panel in x- and y-directions 

in-plane loading parameter Px/^B zero critical flutter speed predicted 
by theory which neglects damping 

static buckling value of 

nondimensional in-plane load parameters; Pjj = Nxa^y^D, Py = Nya^^D 
nondimensional flow speed parameter (see eqs. (24)) 
gas constant 

temperature change from To 
undisturbed air temperature 

undisturbed equilibrium temperature (absolute) of material 

absolute stagnation air temperature 

time 
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t 

u 


nondimensional time , 



undisturbed supersonic airflow velocity 


w lateral displacement 

Xn(?) coefficient functions in a series expansion of (p in trial functions Yn(T?) 

(see eq. (B5)) 

x,y,z Cartesian coordinates (see fig. 1) 

Yn(»?) trial functions in cross -stream direction; free vibration modes of an elastic 

beam 


0! 


coefficient of linear thermal expansion 


degree of thermomechanical couplii^ for plate and beam, respectively; 
1 + V Ea^Tn 


■>'a 


ratio of specific heats for air 


8y2 




flutter determinant corresponding to N -terms of trial function in cross - 
stream direction 


dilatational strain 

C = z/h 
t] = y/b 

6 nondimensional temperature eigenfunction 
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X 


nondimensional eigenvalue 


Xj^ nth eigenvalue, n = 1, 2, . . Re(Xi) > Re(X 2 ) > Re(X 3 ) > . . . 

V Poisson’s ratio 


4 = x/a 


P 


Pa 




<P 

O) 


mass density of panel 

mass density of undisturbed air 

nondimensional thermal diffusion relaxation time for plate, t = WqT* 
nondimensional thermal diffusion relaxation time for rectangular beam of 

pCgh^ / gj \l/2 

width b, thickness h, and length l\ t ,= — ^ 

'’1' UmV 


physical thermal diffusion relaxation time for plate, r* 
nondimensional displacement eigenfunction 
nondimensional eigenfrequency, o) = Im(X) 



/3 \ P^E^^ 
- 2v) k 


Wr 


O) 


Subscript: 


cr 


reference circular frequency, 
circular frequency, w* = 


a>o = 


pha'^y 


denotes value at flutter 


ANALYSIS 

Coupled Equations of Motion 

Consider a fully or partially clamped rectangular panel of length a, width b, and 
thickness h. The panel is situated in a uniform airflow of supersonic velocity U and 
is subjected to the uniformly distributed in-plane compressive forces Nx and Ny 
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(fig. 1). Let T be the temperature change from the constant absolute panel tempera- 
ture Tq of the undisturbed equilibrium state, and let 

ph/2 

Mt = \ Eo:zT(x,y,z,t)dz (1) 

•^-h/2 


be the thermal moment. The temperature change T considered in this analysis is solely 
caused by thermomechanical coupling and, therefore, is small. Hence, the panel material 
constants, such as Young’s modulus E and the coefficient of linear thermal expansion a 
in equation (1), are assumed to be independent of T (but are a function of To). The 
governing equation of small lateral vibrations of the panel is obtained by adding the inertia 


term ph 


at2 


to equation (13.7.1) (with loading 


p = Nxy = Oj of reference 13 (p. 432) as 


follows: 


DV^V^w + N. 


a^w . 

X ^ + -Ny 2" 

3X^ dy^ 


+ ph -^ = F(x,y,t) 

1 - V ^ Qt-S 


( 2 ) 


where 

v2 = A + JL 

8x2 0y2 


(3) 


The aerodynamic force F(x,y,t), according to the two-dimensional quasi-static aerody- 
namic theory, is given by (see ref. 14) 


F(x,y,t) = 


/ 

_ dw 

i , C.o 

~ ax 



(4) 


in which Cg^ is the speed of sound in air, is the air mass density, U is the free - 
stream velocity, and M = U/cg is the Mach number. 
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If the symbols and Ky are used to denote the rotational spring constants 

of the partially clamped edges of the panel at x = 0,a and y = 0,b, respectively, the 
mechanical boundary conditions are 


w = 0 (at X = 0,a; y = 0,b) (5a) 



For simply supported edges, Kx = Ky = 0 in equations (5b) and (5c). For completely 
clamped edges, equations (5b) and (5c) reduce to 9w/9x = 8w/9y = 0. 


Heat -Conduction Equation 

The temperature change T is related to the dilatational strain and the lateral 
displacement w through the coupled heat -conduction equation for plates (eq. (1.12.22) of 
ref. 13, p. 31) given by 


k V^T + 


9z2 / 


3T 

IT " I 


fp 

■ q-Tq — 

- 2v ° at 


( 6 ) 


where k is the coefficient of heat conduction and Cj, is the specific heat at constant 
deformation of the elastic solid. Within the context of the present study (i.e., the elemen- 
tary plate theory of plane stress), the dilatational strain (obtained by using eqs. (8.10.4b), 
(8.10.7), (12.2.1), (12.2.4), and (12.2.6) of ref. 13, pp. 259, 260, 380, and 381) is 


e 


ii " 


2(L^Nt 

(1 - i/)Eh ^ 


1 - 2u 
1 - u 


^2 

zV w 


1 + u 

1 - V 


aT 


1 -2v 
Eh 


(Nx * Ny) 


(7) 


where 

ph/2 

Nrp(x,y,t) = Eq! \ T(x,y,z,t)dz (8) 

^-h/2 

is the in-plane thermal load. Theoretical results for the material damping coefficients 
of an aluminum beam undergoing simple harmonic motion were obtained in reference 11 
for an insulated thermal boundary condition (adiabatic) and an isothermal boundary con- 
dition. These results are compared with experimental measurements in figure 2. As 
can be seen, the experimental results correlate well with the theoretical results obtained 
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o Experimental data (ref. 10) 


-2 



Figure 2.- Correlation of experimental measurements of material damping and theoretical calcu- 
lations of thermoelastic damping coefficient for aluminum beams undergoing simple harmonic 
motion and subject to adiabatic and isothermal boundary conditions. (From ref. 11.) 


for the insulated thermal boundary condition. Hence, the insulated thermal boundary con- 
dition will be assumed here and 


^=0 

dX. 

(at X = 0,a) 

(9a) 

O 

II 

11 

o 

(9b) 

O 

11 

(at z = *|) 

(9c) 


The coupled heat -conduction equation (eq. (6)) can be simplified slightly in the following 
manner: With the aid of equations (7), (8), and (9c), equations (6), (9a), and (9b) can be 
integrated over the thickness (from z = -h/2 to z = h/2) to obtain an auxiliary 
boundary -value problem 


kV^N-p - 



3Eq!^T„\ 8Nt 

1 - 2j/j“9r“ ° 


( 10 ) 


9Np 

Ox 


= 0 


(at X = 0,a) 


(11a) 


9Np 

~w 


= 0 


(at y = 0,b) 


(11b) 
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The solution of equations (10) and (11) is readily obtained as 

N-r = 0 (12) 

Physically, equation (12) indicates an antisymmetric distribution of T over the thick- 
ness. Equation (6) becomes 



g \3T 
1 -2vJ dt 


(1 + 1^)0! V 9t/ 


(13) 


where the nondimensional parameter 


1 H- V 
1 - 


(14) 


has been introduced. This parameter is called "the degree of thermomechanical coupling" 
for a reason which will become apparent later. The coupled equations of motion (eq. (2)) 
and heat conduction (eq. (13)), together with the boundary conditions (eqs. (5) and (9)), con- 
stitute a coupled thermoelastic stability boundary-value problem.. 


Coupled Eigenvalue Equations 
If the nondimensional quantities 


1 = ^ 
^ a 


_ y 


7 ] = 




Vpha^ 


1/2 


(15) 




are introduced into equations (1) to (5), (9), and (13) and the spatial and temporal variables 
are separated by introducing the expressions 

J =?.(£, -,)e« 

^ I 

- y (16) 

To J 

into the resulting equations, the following coupled eigenvalue equations are obtained: 

+ Px -^ + Py(|) + 12(1 + i^)q!To I J dC + q || + (g^x + = 0 


(17) 
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and 


cp(^,r]) = 0 



(at 

4 = 0,1; 

7? = 

0,1) 

(18) 

^ + 12(1 


-.1/2 

'-1/2 

d? T kjj |^= 0 

(at 


{!)) 

(19) 

/a\2 d^cp ^ 
\b/ Qri^ 

12(1 + i^)«To|J 

^ 1 / 2 ® 

^t 

7? = 

0 ) 

(20) 

(if + 



-/3r(h/a) 




(21) 

\a/ 


(1 + 






^=0 




(at 

^ = 

0,1) 

(22a) 

-^=0 

37? 




(at 

7 ? = 

0,1) 

(22b) 

1 Q? 
]] 

O 




(at 

? = 

0,1) 

(22c) 

3^ 









where 


3^2 \b/ 


(23) 


The nondimensional parameters Py, q, g^, and r are associated with the in- 
plane loads Nx and Ny, airflow speed U, aerodynamic damping coefficient and 

the thermal diffusion relaxation time t* , respectively. They are defined by 


Pv = 


Nxa" 


N„a^ 


P =_I 

y D 


q = 


PaU^ 


d)|m^ - 1 


Sa 


D 


T = COqT- 


(24) 




where 


. Cl + — 

\ 1 - 2vJ k 
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and where 


D 


= 


pha^ 


1/2 


( 25 ) 


is the reference circular frequency. The thermal diffusion relaxation time t* is asso- 
ciated with the time required for any temperature perturbation to return to the undisturbed 
equilibrium temperature as a result of thermal diffusion. Equations (17) to (22) govern 
the associated eigenvalue problem of the coupled thermoelastic -aeroelastic boundary- 
value problem where X is the eigenvalue. 


Reduction of Coupled Eigenvalue Equations to 
Uncoupled Eigenvalue Equations 

For a thin panel, the thickness -length ratio h/a is much smaller than unity. 

2 2 

Accordingly, the in-plane heat flow term V 0 in the heat -conduction equation (21) 
can be neglected, since its effect on the temperature eigenfunction 6 is much smaller 
than that of the heat flow term in the thickness direction. If the in -plane heat 

flow term in equation (21) is neglected, the edge thermal boundary conditions (eqs. (22a) 
and (22b)) cannot be prescribed arbitrarily; that is, they must be ignored. The simpli- 
fied heat -conduction equation, together with the remaining thermal boundary conditions 
(eq. (22c)), is then solved for 9 in terms of the unknown displacement eigenfunction cp, 
and the solution is 


ea,v,0 


g(h/a)V^ffl 


«To(l + u)(^ 


1 + 



sinh\/^g 


cosh 


\/^l 


(26) 


With 9 known as a function of <p, equations (17), (19), and (20) can be integrated to give 
the followii^ uncoupled governing equations and boundary conditions for <p: 



a cp 

9772 




)<p 


(27) 


= 0 

(1 * g)% 



= 0 




877 


(at I = 0,1) (28a) 


(at 77 = 0,1) (28b) 


12 



where 



The parameter g is the thermomechanical coupling function for thin plates since it 
represents the effect of thermomechanical coupling on the dynamic system. The real 
part of g(Xr) (where X is usually complex) represents a stiffening effect, whereas the 
imaginary part of g(Xr) represents a damping effect. Since j3 appears as a factor in 
the expression for g, it is appropriate to refer to this parameter as the degree of thermo- 
mechanical coupling. 

Equation (27) reveals that material damping (thermomechanical coupling function g) 
modifies only the bending term of the governing equation. Similarly, flutter equations 
derived on the basis of viscoelastic theory utilizing stress -strain relations for a Kelvin- 
Voigt body indicate that linear hysteretic structural damping should modify only those 
terms of the equations associated with bendir^ (ref. 6). As pointed out in reference 6, 
many authors (refs. 15 to 20) have introduced a hysteretic structural damping coefficient 
that, in effect, is used to modify both the bending and membrane loading terms of the flut- 
ter equation. It can be shown by writing the equation of energy balance that these authors 
have unknowingly introduced negative membrane damping into their analyses. A com- 
parison of the equations derived herein on the basis of thermoelasticity with those pre- 
sented in reference 6 based on linear hysteretic damping reveals that the latter develop- 
ment neglects the stiffening effect associated with the real part of g(Xr) and also com- 
pletely neglects the effects of g in the rotational boundary conditions. 


Solution Technique and Flutter Criterion 

The uncoupled eigenvalue equations (27) and (28) cannot, in general, be solved 
exactly if ky 0. Therefore, one must resort to an approximate technique such as 
Galerkin's method (e.g., see ref. 21) for a solution. In the present investigation the solu- 
tion method uses a generalized Galerkin's procedure (wherein only geometric boundary 
conditions must be satisfied) in the cross -stream (rj) direction to reduce the governing 
partial differential equation to a system of ordinary differential equations in the stream- 
wise (^) direction that approximately governs the problem. The latter equations are 
solved exactly. This solution method is quite general and can be applied to various panel 
flutter problems including, for example, the effects of arbitrary airflow direction (refs. 21 
and 22). This method is more suitable than the procedure where the Galerkin method is 
used in both the cross -stream and streamwise directions because, with the present method, 
approximations are involved in one direction only. For a given accuracy, the present 
method usually leads to a smaller flutter determinant. The sizes of the determinants for 
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the present and Galerkin method are (4N) x (4N) and (N X J) X (N x J), respectively, 
where N and J are the numbers of terms in the trial functions in the cross -stream 
and streamwise directions, respectively. Details of the solution are given in appendix B. 

NUMERICAL RESULTS AND DISCUSSION 

Thermoelastic Damping 

The effect of thermoelastic damping on the coupled thermoelastic -aeroelastic sys- 
tem is represented by the imaginary part of the thermo mechanical function g(Xr) (cf. 
eq. (29)) and is given by 

g„[Re(X)r,a;T,j3,z^] = ^ — Im 1 - ^ + — tanh ^ (w = Im(X)) (30) 

1 - 2v 

The quantity grp is known as the "thermoelastic damping coefficient." The term in 
brackets is a function of cut (cut = cu*t* where cu* is the circular frequency, 
cu* = the definitions of /3 (eq. (14)) and t* (eq. (24)), g,p at flutter is 

seen to be a function of the panel material properties E, v, a, pCg, and k, the abso- 
lute equilibrium temperature Tq, the panel thickness h, and the panel flutter frequency 
cu*j,. The effect of panel length a and width b on grp is implicit through their effect 
on the flutter frequency; grp is not an explicit function of a and b. 

St* 

The variation of the ratio — (at flutter) with cut, as obtained from 

equation (30), is shown in figure 3. As <^r increases from 0 to the ratio increases 
1.0 
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from 0 rapidly and monotonically to a maximum value of 0.494 at wt = 9.9; thereafter, 
it decreases monotonically to 0 at cot = +«>. The material properties for 2024 -T3 alu- 
minum alloy at room temperature are as follows: 


p 2705 kg/m3 (5.25 slugs/ft3) 

E 72.4 GN/m2 (10.5 X lO^ psi) 

a 23 X 10-6 per K (12.8 X 10"6 per op) 

Ce 0.84 kJ/kg-K (0.20 Btu/lb-OF) 

k 0.126 kW/m-K (0.243 Btu-in./ft2-s-°F) 

V 0.33 


With these properties and under a room temperature Tq = 300 K (540° R), the degree of 
thermomechanical coupling jS is calculated from equation (14) to be /3 = 0.01, and from 
figure 3 the maximum value of grp is found to be 0.0048. 

St* 

The trends exhibited by the variation of the ratio — ± with frequency 

7 l~- 2 ^ 

for thin plates (fig. 3) are seen to be the same as that for beams with rectangular cross 
sections (see fig. 2). However, it should be noted that the experiments of Bennewitz and 
Rotger (ref. 23) revealed that for very low frequency oscillations the value of material 
damping for wires may continue to decrease monotonically with decreasing frequency or 
begin to increase depending on the material; material damping for steel, brass, and glass 
decreased but results for silver and aluminum showed increases. Since that part of 
material damping due to plastic flow, which is neglected in the present analysis, is larger 
the lower the frequency, the increases in material damping at very low frequencies were 
attributed to plastic flow. This seems plausible since very accurate measurements have 
revealed the presence of permanent strains for aluminum specimens for very small 
stresses, whereas, for mild steel specimens, such permanent strains (if any) may be 
negligible (ref. 24). The experimental results of reference 23 reveal that the peak value 
of the material damping coefficient for aluminum wires is 2.6 x 10“3, which occurs at 
COT = 83 Hz. These results agree extremely well with the theoretical value of the damping 
coefficient of 2.3 x 10“3 at cot = 84 Hz for thermoelastic damping (ref. 25). 

The good correlation between theoretical and experimental damping results for 
aluminum beams (fig. 2) and for wires (refs. 23 and 25) suggests that material damping 
for thin aluminum panels, at least in the rai^e of frequency for which the effects of plas- 
tic flow are negligible, consists mainly of thermoelastic damping. Hence, material 
dampii^ coefficients for aluminum panels are for practical purposes given by equation (30) 
if COT is not vanishingly small. 

Measurements of structural damping for panels (e.g., ref. 26) are usually consider- 
ably greater than the thermoelastic material damping coefficients given by equation (30). 
However, these measurements include not only material damping but also other forms of 
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damping such as frictional damping at the panel supports and, in the case of built-up 
panels, the contribution of other forms of damping may be considerable. 

In concluding the discussion of thermoelastic material damping, it is suggested that 
there is a need for additional experimental effort on structural damping for beams or 
panels, particularly for small values of wt (wt < 10). 

Flutter Boundaries 

The solution procedure outlined in appendix B for finding an approximate flutter 
speed parameter q^r has been programed for a CDC 6600 electronic computer. The 
functions Yn(^) in the trial function cpU,r]) given by expression (B5) are chosen to be 
those of free vibration modes of the elastic beam satisfying the boundary conditions (28b) 
with g = 0 (see appendix C). Numerical results for flutter speeds in terms of 

ql/3 and flutter frequencies Wgr are presented in figures 4 to 7 for flat panels of 
2024-T3 aluminum alloy with thickness h = 0.13 cm (0.053 in.), thickness -length ratio 
h/a = 0.00204, and aspect ratios of b/a = 0.346 and 0.707. The absolute undisturbed equi- 
librium temperature Tq is assumed to be that of room temperature Tq = 300 K (540° R) 
and the undisturbed airflow condition is assumed to be at a Mach number of 3 and an abso- 
lute stagnation temperature Ts = 428 K (770° R). These panel configurations (with 
b/a = 0.346) and flow conditions correspond to a test case given in reference 27. For 
these conditions and the material properties given previously for 2024 -T3 aluminum alloy, 
the nondimensional thermal diffusion relaxation time t is calculated to be 0.157, the 
aerodynamic damping coefficient g^ is calculated to be 0.00405q (see appendix D), and 
the degree of thermomechanical coupling ^ is calculated to be 0.01. 

All flutter results are based on the two-term solution (N = 2 in eq. (B5)) by using 
the symmetric first and third free vibration modes of the elastic beam. Since the second 
mode is antisymmetric, inclusion of this term in the analysis has no influence on the 
flutter results given herein. While differences between the one- and two-term solution 
for the critical speed parameter q^/^ are observed to be negligibly small for the com- 
pletely clamped panel, the difference becomes somewhat significant for partially clamped 
panels. The maximum difference for ql/3 corresponding to the case of curve 2 in fig- 
ure 6 is 2 percent (6 percent for q itself). 

Figure 4(a) shows the theoretical flutter boundaries for the completely clamped 
panel for undamped and damped cases. The cases for which no damping or aerodynamic 
damping alone is present were studied in references 14 and 28. For the undamped case 
(curve 1), the flutter speed q^r becomes zero at Px/Pb = 0.94 s P, where Pg is the 
static buckling load for q = 0. Since the aerodynamic damping coefficient g^ is pro- 
portional to q, the presence of aerod 3 mamic damping (cf. curve 2) has little influence on 
the flutter boundary for no material damping (|3 = 0) in the vicinity of P, although it does 
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remove the physically untenable condition of flutter at q = 0. The effect of thermoelastic 
damping (curve 3) is seen to be of considerable importance in the vicinity of P, although 
it has little influence on the flutter boundaries for smaller values of Px/Pb and no 
influence at the intersection point of the flutter and buckling boundaries where o) =5 0 
and thus g>j< = 0. 

The flutter frequencies corresponding to curve 3 of figure 4(a) are plotted in fig- 
ure 4(b) as a function of the in-plane load parameter Px/Pg. The corresponding ther- 
moelastic dampir^ coefficient grp is given in figure 5 as a function of the flutter fre- 
quency o)cr* As the in-plane load parameter Px/Pb increases, flutter frequency Wcr 
decreases (fig. 4(b)), but the damping coefficient grp increases initially with decreasing 
frequency (fig. 5) and then decreases to zero at = 0, which corresponds to the inter- 
section of the flutter and buckling boundaries. The value of grp reaches its maximum 
value in the vicinity of Px/^b ~ ^ (where q^j. = 0 when there is no damping). 



Figure 5.- Variation of thermoelastic (jamping coefficient with flutter frequency 
for the compfete/y clamped panel considered in figure 4. 


Figure 6 shows the flutter boundaries for a partially clamped aluminum panel with 
kx = ky = 30 for the panel configuration and panel equilibrium temperature considered 
in figure 4 for the completely clamped panel. Results are presented for no damping and 
for both aerodynamic and thermoelastic damping. (See also ref. 29 for results of 
undamped cases corresponding to various values of kx and ky.) The zero flutter speed 
observed for the completely clamped panel for no damping no longer appears in the pres- 
ent case, and the effect of damping on panel flutter is seen to become less important than 
for the fully clamped case. However, if the aspect ratio for the present panel is allowed 
to vary, it is again possible to find a case for which q^r = 0 absence of all kinds 

of damping. 
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Figure 6.- Flutter boundaries for a partially clamped aluminum panel, b/a = 0.346; h/a = 0.00204; 
h = 0.13 cm (0.053 in.); Py/Px = 1; t = 0.157; kx = ky = 30. 

Figure 7 shows the flutter boundaries for a simply supported panel^ with the same 
airflow condition, panel equilibrium temperature, and panel configuration as considered 
previously except that the aspect ratio b/a is now purposely chosen to be 0.707 so that 
for no damping the flutter speed is zero at Px/Pb* effect of aerodynamic and 

thermoelastic damping is seen to be small and destabilizing, ^ in particular, damping has 
no effect at Px = Pb since the flutter frequency is zero for this condition and both aero- 
dynamic and thermoelastic damping are frequency dependent and vanish as the frequency 
approaches zero. Hence, it is apparent that inclusion of only aerodynamic and thermo - 
elastic damping in the theoretical flutter analysis cannot completely remove the para- 
doxical phenomenon of zero flutter speed predicted by the theory with no damping for all 
ranges of edge restraint and aspect ratio. 

In order to remove the paradoxical phenomenon, one must include frequency- 
independent hysteretic material damping and support damping. An example of the exis- 
tence of frequency-independent support damping is the in-plane damping caused by the 
in-plane kinetic dry friction force between the panel and support contacting surfaces (see 
refs. 32 and 33). An alternative approach is to account for possible initial imperfections 



^Equations (27) and (28) can be solved exactly for simply supported panels. (See 
refs. 30 and 31 for exact solution with g^ = 0.) 

^As shown in reference 11, small amounts of thermoelastic damping may have a 
destabilizing effect on a nonconservative system. 
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Figure 7.- Flutter boundaries for a simply supported aluminum panel, 
b/a = 0.707; h/a = 0.00204; h = 0.13 cm (0.053 in.); Py/Px = 0; 

T = 0.157; P = 0.01; = 0.00405q. ^ 

in the flutter analysis (see ref. 8). This approach, as well as the case of in -plane sup- 
port damping, requires a nonlinear analysis, which is beyond the scope of the present 
study. 


CONCLUDING REMARKS 

Thermoelastic material damping and its effect on the flutter of thin rectangular 
panels subjected to in-plane stresses in a supersonic airflow has been studied within the 
context of linear thermoelasticity theory. A closed-form expression for the coefficient 
of thermoelastic damping was obtained as a function of frequency, panel thickness, panel 
temperature, and material properties. The solution procedure for solving the flutter 
boundary -value problem consisted of application of a generalized Galerkin technique in 
the cross -stream direction. This reduced the governing partial differential equations to 
a set of uncoupled ordinary differential equations in the streamwise direction and these 
equations were then solved exactly. 

Numerical results for the thermoelastic material damping coefficient and flutter 
boundaries for simply supported and partially or fully clamped panels were presented. 

It was shown that the importance of the effects of thermoelastic damping on flutter bound- 
aries depends on the panel configuration and boundary conditions. Calculations were pre- 
sented for a panel with a zero flutter speed (for no damping) which occurred at a frequency 
for which the thermoelastic damping coefficient was near its peak value. For this condi- 
tion the effect of thermoelastic damping was shown to be quite significant. However, 
thermoelastic damping had no effect at the point where the flutter and buckling boundaries 
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intersect (for which the flutter frequency is zero) and may even have a destabilizii^ 
effect for certain conditions. The cases studied indicated that inclusion of thermoelas- 
tic dampii^, in general, may remove the anomalous behavior of zero flutter speeds pre- 
dicated by flutter theory for no damping. However, this is not true when flutter speed 
occurs at a value of in-plane stress required for bucklir^. In order to remove such an 
anomalous behavior completely, frequency-independent structural dampii^ must be 
included or effects of initial imperfections must be considered. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., August 27, 1971. 
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APPENDIX A 


CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 

Factors required for converting the units used herein to the International System 
of Units (SI) are given in the following table: 


Physical quantity 

U.S. Customary 
Units 

Conversion 

factor 

SI Unit 



(*) 


Density 

slugs/ft^ 

515.379 

kg/m^ 

Force 

Ibf 

4.448 

N 

Gas constant 

ft2/s2_OF 

0.1672 

m2/s2-K 

Heat energy 

Btu 

1054 

J 

Length 

(r 

0.0254 

]” 

\ft 

0.3048 

j 

Mass 

fslugs 

14.59 

jkg 


|lbm 

0.45359 

Specific heat 

Btu/lbm-°F 

4184 

J/kg-K 

Stress and pressure 

psi = Ibf/in^ 

6895 

N/m2 

Temperature 

fOF 

l°R 

(5/9)(F + 459.67) 
5/9 

y 

Thermal conductivity 

Btu-in./ft^ -s -Of 

518.87 

W/m-K 


*Multiply value given in U.S. Customary Unit by conversion factor to obtain equiv- 


alent value in SI Unit. 

refixes to indicate multiples of SI units are as follows: 


Prefix 

Multiple 

giga (G) 

10^ 

kilo (k) 

10^ 

centi (c) 

10-2 
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APPENDIX B 


SOLUTION OF THE UNCOUPLED EIGENVALUE EQUATIONS 


The generalized Galerkin method (ref. 34) is employed to reduce the governing 
partial differential equations to a set of ordinary differential equations that approximately 
govern the problem. Multiplication of equations (27) and (28a) by d<p, the first variation 
of and integration of the resultii^ equations over 77 from 77 = 0 to 1, using the inte- 
gration by parts technique and the boundary conditions (28b), yield 




where 


+ X(X + g^) (p5(p d77 + h(4) = 0 

r (77(0,77) 6 <;p d77 = C <77(1,77)6(77 d77 = 0 

-•o '^0 

(1 + g) r T kx 1^ 6^ d77 = 0 


h(^) 


^ J'k,[^{,0) Mi * ^U,l) (ky . •») 

(ky = 00 , completely clamped edges^ 


(Bl) 

(B2) 

(B3) 

(B4) 


An approximate solution of equation (27) is sought in the form 
N 

<pU,V) = ^ X„(?)Yn(t7) (N = Integer g 1) (B5) 

n=l 

where Yn(T?) is a set of (given) known relatively complete functions (see ref. 35) and 
^( 4 , 77 ), as stipulated by the principle of virtual work, must satisfy all geometrical boundary 
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APPENDIX B — Continued 


conditions, that is, zero displacement boundary conditions in equation (28b) if ky °o 
and both zero displacement and zero slope boundary conditions in equation (28b) if ky = <» 
(completely clamped edge conditions). Since Yn(?7) is a set of known (given) function, 
there is no variation on Yniv) (that is, SYniv) = 0 for n = 1, 2, . . ., N) and the varia- 
tion of the function cp is 

N 

^(pU,'n) = Yj 5Xn(^) (B6) 

n=l 


Substitution of equations (B5) and (B6) into equations (Bl) to (B3), since the variation 6Xn 
is arbitrary, leads to 


N r 

Y S)amnXn (^) - 2(1 + bj^n " ^x^mn ^n(^) + *l^mn^n(^) 

n=l 


(1 + Cmn - Pyt>mn(|) + ^)^mn + (|) 


Xn(^) > = 0 (B7) 


N 


n=l 

N 


N 

I 

n=l 


J a„nX„(0)= ^ a„„X„(l) = 0 


+ g)x”( 0 ) -kxX;( 0 )] =0 


(m = 1, 2, . . ., N) (B8) 


(B9) 


n=l 


where 


N 

[(1 + g)x;’(i) + kxx;(i[| = 0 


n=l 


amn = \ Yn(??)dr? 

-'O 

p 1 

Vn = J^ Y^(?7) Yn(r?)dr? 

c ^ 

Y^(7) Yn(?7)d?7 


^mn 


*mn 


fky[Ym(0) Y» + Y^d) Y;(d] 


(BIO) 


(ky * °o) 
(ky = -) 
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APPENDIX B - Continued 

The general solution of the system of N-homogeneous linear differential equations (B7) 
is of the form 


Xn(4) = J Anje^i^ 
i=i 

Substitution of equation (Bll) into equation (B7) yields 


(n = 1, 2, . . N) (Bll) 


N 

^ Ajjj = 0 (m = 1,2,. . ., N; n = 1, 2, , . ,, 4N) 

i=i 

where 


Qmn(^i) - (1 + g) 2-mn^j 


2(1 +g)Vn^' 




l*x^mn 


^ (1 + g) Cmn ^ l*yl^mn 


3 


+ X 


(^A + 




hmn^ 


and S2j are the 4N-roots of the determinantal equation 


(B12) 


(B13) 


det Qmn(^) = ® = 1, 2, . . ., N; j = 1, 2, . . ., 4N) (B14) 

If Aj^^ is denoted by An, then equation (B14) can be solved for Ajn (j = 2, 3, . . ,, N) 
in terms of An to give 

Ajjj = fj(^^)An (j~l>2,. . ., N) 

fl(^i) = 1 

Thus, expression (Bll) becomes 
4N 

Xj(^) = ^ fj(l2n)Ane“>^^ (i = 1, 2, . . ., N) (B16) 

1^1 

Finally, expression (B16), together with the 4N -boundary conditions (B8) to (BIO), pro- 
vides the following 4N -linear homogeneous algebraic equations in An: 

amjfj(^n)An = 0 (B17) 


amilj(«n)e^*"An = 0 (m = 1, 2, . . ., N) (B18) 


4N N 

n 

n=l j=l 
4N N 

zz 

n=l 3=1 
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APPENDIX B — Continued 


4N N 

^ ^ ajnjd +g)fj(On)fin(«n -kx)An = 0 (B19) 

n=l j=l 

4N N 

^ ^ + kx^e = 0 (B20) 

n=l j^l 

For a nontrivial solution, equations (BIT) to (B20) yield the characteristic equation 

AN[^>q.gA’^x>Py»g(/3,!^,XT),b/^ = det H„in = 0 (B21) 

where Hmn is the 4N x 4N coefficient matrix of An in equations (BIT) to (B20). 

Since equation (B21), in general, is a transcendental equation in A, it yields infinitely 
many eigenvalues (roots) for X, 

>-n = >-n(<l>gA>Px»Py>/3,T,b/a) (n = 1, 2, . . .) (B22) 

as functions of q, g^, Px, Py, t, and b/a. In the presence of damping these 
eigenvalues are, in general, complex. Let Xn (n = 1, 2, . . .) be ordered in such a way 
that 

Re(Xi) £ Re(X2) ^ Re(X3) . . . (B23) 

where Re^Xnj stands for "real part of Xn." The undisturbed equilibrium state is then 
said to be (see ref. 36) asymptotically stable if all the real parts of Xn are negative 
^i.e., Re(Xi) < O), dynamically unstable (flutter) if one or more real parts of Xn are 
positive and the corresponding imaginary parts of Xn (frequencies) are nonvanishing, 
and statically unstable (divergence or buckling) if one or more real parts of Xn are 
positive and at least one of the corresponding imaginary parts of Xn vanishes. 

For a given set of values of g^, jS, t, b/a, v, Py/Px, and Px smaller than 
its static buckling value Pg, equation (B21) may be written as 

AN(^,q) = 0 , (B24) 

For the numerical solution of this equation, the followii^ method was used: The flow 
speed parameter q is gradually increased from 0 and the eigenvalue X corresponding 
to the lowest frequency (lm(X)) among all eigenvalues Xn is calculated from the charac- 
teristic equation (eq. (B22)). The eigenvalue so obtained must have a negative real part 
in the region of stability which, for the present problems, corresponds to q < q^,j.. The 
flutter speed q^^ is reached when q is increased to the value for which the real part 
of the eigenvalue vanishes and further increases in q result in positive real parts. If 
the instability is associated with the second lowest frequency, the real part of the eigen- 
value associated with the lowest frequency becomes discontinuous, changes sign from 
negative to positive abruptly, and does not vanish for any value of q. For this case, one 
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APPENDIX B - Concluded 


would calculate the eigenvalue associated with the second lowest frequency rather than 
the first and find q^j. so that the real part of this eigenvalue vanishes. However, for 
all calculations presented herein for panels, the instability was associated with the eigen- 
value corresponding to the lowest frequency. 


27 


APPENDIX C 


FREE VIBRATION MODES OF AN ELASTIC BEAM 


where C is an arbitrary constant, 


(n = 1, 2, . . .) 


6n = 


sin a 

cosh f2„ - cos f2n + ; 

“ ky b 

sinh fin - sin fi- 


~ 0X11 

and fin = Im(^n) is the nth root of the characteristic equation 

( 2fin a\ 

cosh fin + cos fin - 5n sinh S^n ” ®n + “i — - sin fin 

V % 7 _ 

^ / 2 fi« «\ 

- ky Sinh fin + sin fin - 6n cosh fin + 1 6n + ^n 


= 0 


^fij^ < fi2 < fi0 


(Cl) 


The governing equations of free vibration of an elastic beam with a nondimensional 
length of unity and rotational spring constant ky are obtained from equations (27) 
and (28) by setting cp = Y(? 7 ) and Px = Py = g = q = = 0 as follows: 

Y""(?7) + x2y(?7) = 0*^ 

Y(0) = Y(1) = 0 

Y”(0) -|kyY(0) = 0 

Y”(l) +|kyY(l) = 0 

The solution of equations (Cl) for ky # 0 is easily found to be 

r /2fi„ a \ 

Yn(’l) = C cosh fin^ - cos fi^?? - 6n sinh ^ + 6jysin fin?7 


(C2) 


(C3) 


(C4) 
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APPENDIX D 


AERODYNAMIC DAMPING COEFFICIENT 


From the definitions of q and in equations (24), the aerodynamic damping 
coefficient is given as follows: 


SA = 




CaM^a 



where the relation M = U/ca has been used. Since for air 


(Dl) 



where the ratio of specific heats for air equals 1.4, the gas constant for air R 
equals 287 m2/s2-K (1717 ft2/s2-OF), and Tg and Ta are the absolute stagnation tem- 
perature and the undisturbed air temperature, respectively, equation (Dl) reduces to 
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